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We report a cluster of results regarding the difficulty of finding approximate ground states to 
typical instances of the quantum satisfiability problem fc-QSAT on large random graphs. As an 
approximation strategy, we optimize the solution space over 'classical' product states, which in 
turn introduces a novel autonomous classical optimization problem, PSAT, over a space of contin- 
uous degrees of freedom rather than discrete bits. Our central results are: (i) The derivation of 
a set of bounds and approximations in various limits of the problem, several of which we believe 
may be amenable to a rigorous treatment, (ii) A demonstration that an approximation based on 
a greedy algorithm borrowed from the study of frustrated magnetism performs well over a wide 
range in parameter space, and its performance reflects structure of the solution space of random 
fc-QSAT. Simulated annealing exhibits metastability in similar 'hard' regions of parameter space, 
(iii) A generalization of belief propagation algorithms introduced for classical problems to the case 
of continuous spins. This yields both approximate solutions, as well as insights into the free energy 
'landscape' of the approximation problem, including a so-called dynamical transition near the satis- 
fiability threshold. Taken together, these results allow us to elucidate the phase diagram of random 
fc-QSAT in a two-dimensional energy-density-clause-density space. 
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I. INTRODUCTION 

Ever since the realization that quantum mechanics 
could be used to solve computationally difficult problems, 
there has been an intense effort to understand quantum 
computers and quantum algorithms. Naturally, much 
of the algorithmic research front has been focused on 
those that give a substantial speed up over their classical 
counterparts. Complexity theory suggests, however, that 
there are many natural optimization problems that are 
hard even for a quantum computer to solve. These in- 
clude both classical NP-complete problems such as satis- 
fiability (SAT) and intrinsically quantum QMA-complete 
problems such as quantum satisfiability (QSAT) which 
require a quantum computer even to represent the rele- 
vant state space. These hardness results suggest the need 
for approximate strategies where one settles for less-than- 
globally-optimal solutions. Unfortunately, this too may 
be hard as evidenced by the classical PCP theorem [HE]- 

The PCP theorem, in this setting, implies that there 
are classical energy functions for which it is NP-complete 
not only to find the ground state but even to find highly 
excited states with an energy density below some thresh- 
old. This result has astonishing dynamical consequences. 
Under the assumption that nature cannot solve NP- 
complete problems, there are no dynamics which can 
relax these systems to low temperature. This consti- 
tutes an abstract argument for the existence of glassiness 
which, in an amusing parallel with the usual statistical 
analysis of spin glass theory, makes no reference to the 
specific dynamics of the system. 

As we expect that quantum computers find NP- 
complete problems just as troubling as classical comput- 
ers, the classical problems that are hard to approximate 



remain so even exploiting quantum dynamics. However, 
this leaves open the question of how hard intrinsically 
quantum (QMA-complete) problems (HHi] might be to 
approximate, i.e. of the precise statement of a quantum 
analog of the PCP theorem. One possibility derives from 
the usual intuition of quantum statistical mechanics that 
a finite energy density corresponds to a finite tempera- 
ture which renders quantum mechanics essentially clas- 
sical. If this intuition governs then we would expect to 
be able to efficiently represent approximate finite energy 
density states for QMA-complete problems on a classi- 
cal computer and thus the approximation problem would 
"only" be in NP and thus NP-complete. The other, more 
interesting, possibility is that the above is not true and 
that there exist quantum Hamiltonians whose low energy 
density states are both hard to find and impossible to rep- 
resent with an efficient classical description. Here even 
the approximation problem would be in QMA. Of course 
these possibilities are not exclusive — possibly there is a 
"phase transition" as a function of the energy density 
where the problem goes from NP to QMA as the energy 
density is lowered. 

In this paper we take an approach to this question 
which fits naturally within the quantum statistical me- 
chanics of random systems wherein, unlike in the Com- 
puter Science setting, we do not attempt to make state- 
ments about all possible instances of a problem but in- 
stead pick a natural measure on the space of problems 
and try to investigate the properties of instances which 
are typical with respect to it. Specifically, we consider 
the approximation problem for the low energy states of a 
canonical QMA-complete problem: /c-body quantum sat- 
isfiability (fc-QSAT). As in previous work [HHS], we con- 
sider a uniform random ensemble of fc-QSAT instances 
with M clauses, N qubits and clause density a = M/N 
held fixed in the thermodynamic limit. 

The current state of knowledge and our general pro- 
gram with regard to this ensemble is summarized in Fig.fl] 
where we consider a "phase diagram" in a plane labeled 
by a and the energy per clause. It is known that typi- 
cal instances are satisfiable (i.e. have strictly zero energy 
ground states) by product states at low density, a < aps , 
satisfiable by entangled states for intermediate densities 
aps < a < as, and strictly unsatisfiable for high den- 
sities a > as- For fc = 3, of maximum interest in this 
paper, Ups =0.91... and a^ = 1.0 ±0.06 [71l2].[Tn]. The 
questions we would like to answer are denoted by a set 
of conjectured "phases" in the Fig. 1. Specifically, these 
are first to establish the boundaries of the regions where 

• it is possible to use product states to achieve the 
required energy densities 

• it is possible to use non-product classical states to 
achieve the required energy densities 

• it is essential to use quantum states to achieve the 
required energy densities 
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FIG. 1. Heuristic 'hardness' phase diagram from approximat- 
ing random QSAT. PS regions have product state witnesses, 
C have classical non-product and Q have only quantum wit- 
nesses. "Easy" and "Hard" refer to conjectural algorithm- 
dependent hardness transitions for finding such witnesses. 



• it is not possible to achieve the specified energy 
densities 

By "classical states" we mean states for which energies 
can be evaluated in a time poly(A^) by a classical com- 
puter. Additionally, there is the question of how long 
it takes to find such states. In principle for each cat- 
egory of states there can be algorithm dependent easy 
and hard regions. It is also a challenge to determine 
these boundaries. Needless to say, it is not known a 
priori whether all of these regions actually exist (apart 
from the product state and inaccessible regions). We 
do know though that the random ensemble exhibits a fi- 
nite surface/volume ratio for subgraphs of its core; this 
'expansion' property makes it resistant to simple divide- 
and-conquer approximation strategies that apply to finite 
dimensional graphs and thus leaves open the possibility 
of hard-to-approximate phases. 

In this paper we take a step towards answering the 
above questions by attempting to find the boundaries for 
the product state regions. Working over the variational 
class of product states for QSAT defines a classical, con- 
tinuous spin, energy functional on the same interaction 
graph. We dub this the PS AT model (see Eq. ^). The 
analysis of this model is the central technical goal of this 
paper. By introducing a classical Gibbs ensemble with 
temperature T = 1//3 onto PSAT, we can use statisti- 
cal mechanical techniques adapted from cavity analysis 
[TTl [T^ to probe the structure of the finite energy den- 
sity landscape of this classical model. We note that cav- 
ity analysis of the associated classical model is similar in 
spirit to the cavity analysis of the coherent state represen- 
tation of the AKLT spin glass [131 El j but the quantum- 
classical mapping involved here is approximate. We will 
extract complementary information regarding the low en- 
ergy structure by direct numerical optimization of finite- 



size instances using both a local greedy quench algorithm 
and simulated annealing. 

Using these approaches, our aims are two-fold. First, 
we extract the ground state energy density of PSAT as a 
function of density a, as this defines the energy density 
below which QSAT has no product state approximation. 
Second, we accumulate evidence for the existence of a 
classically 'hard' phase at energy densities approaching 
the ground state boundary from above. This phase ex- 
hibits metastability in the quench dynamics at finite size 
and a 'dynamical replica symmetry breaking' instability 
in our cavity analysis |11| . In passing, we will provide a 
characterization of the low and high temperature ther- 
modynamics of PSAT. 

Finally we note that our work here is complementary 
to more rigorous partial results on the quantum PCP 
problem. The development of a quantum gap amplifica- 
tion procedure [15] provided a quantum analogue of the 
central step in Dinur's proof of the classical result [IB] . 
However, the quantum no-cloning theorem has stymied 
attempts to complete a quantum PCP proof along those 
lines. Rather, the field has turned to constructing clas- 
sical approximation strategies and classical approximate 
witnesses to restricted classes of Hamiltonians in an at- 
tempt to constrain possible quantum PCPs. For exam- 
ple, a now classic result is that commuting 2-local Hamil- 
tonians have only short-range entangled ground states 
[T7] which possess efficient classical representations and 
no topological order. Several workers have extended 
these considerations to finite energy density witnesses to 
some classes of fc > 2-local commuting Hamiltonians [TSl - 
I21j under various restrictions. Of particular relevance to 
us is the result that commuting Hamiltonians wherein a 
small fraction of the terms can be dropped to render the 
interaction complex 1-localizable possess classically rep- 
resentable states of low energy density [20j [21] . We note 
that the random ensemble we study in this work is in- 
deed 1-localizable but not commuting. Ground states 
constructed for commuting Hamiltonians may be per- 
turbatively extended to finite energy density states for 
almost-commuting Hamiltonians [55]. Finally, for gen- 
eral (non-commuting) quantum Hamiltonians, the rigor- 
ous results have relied on constructing guarantees that 
product states [23j or coarse grained product states (in 
the 2-local case) [21] exist which provide good approxi- 
mations to the ground state energy under certain geomet- 
ric constraints on the interaction graph. Of these results, 
we note that low degree expanders as we study in this pa- 
per are precisely in the class that evades existing efficient 
approximation guarantees. 

In the balance of the paper we do the following. We 
describe in the next section the random QSAT model, 
its relationship to the classical PSAT model and various 
pieces of useful terminology. Section [TTTj discusses various 
limiting expressions for the PSAT thermodynamic quan- 
tities. Some details of these calculation can be found in 
the appendices. These analytic results serve as a useful 
test of the numerics. Next, we present results from our 



greedy local quench, simulated annealing and exact di- 
agonalization algorithms in Section IV In Section [Vj we 
analyze the product state phase with the cavity method 
where we find some evidence for a dynamical instabil- 
ity. Finally, we give concluding remarks in Section VI 



As an aid to the reader, we collect several results on the 
Haar measure on the space CP"~^ in Appendix JBJ and 
hypercores in Appendix [C] 



II. MODEL 



The QSAT Hamiltonian on N qubits is given by 



H 



En- 



(1) 



where H™ is a /c-local rank 1 projector associated with 
the (hyper-)edge m — (7711,7712, •• • ,mk) of an interac- 
tion graph G. H™ projects the state of k qubits onto one 
'direction' in their 2*''-dimensional joint Hilbert space, ex- 
acting on energy cost given by the size of the projection. 
Such a Hamiltonian has a zero energy state |^) if and 
only if 1^) is simultaneously annihilated by all of the pro- 
jectors H™. In this case, we say that H is satisfiable and 
that each of the clauses H™ is satisfied by the state |^). 
We note that deciding whether a given QSAT Hamilto- 
nian H is satisfiable is QMAi-complete for A: > 3 [51l25|. 
subject to the technical restriction that the ground state 
energy of H , if non-zero, is larger than a polynomially 
small promise gap S = l/poly(N). As we expect that 
the energetics of QSAT in the ensembles we consider are 
extensive, the promise gap should be satisfied with high 
probability and we will not worry about it further. 

The classical PSAT Hamiltonian on N vector spins is 
given by the variational energy of the QSAT Hamiltonian 
H over the class of qubit product states: 



Hi{h,})=Y,{{f^^}^e^m\Tl"'\{n,} 
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dm) 



(2) 
(3) 



Here, \fii) is the spin-1/2 coherent state pointing in the 
direction of the unit vector fii at site i and |0™) is the 
fc-local state onto which the projector H™ projects. In 
complex coordinates, the Hamiltonian takes on the k- 
quadratic form 



H = 



/ . IV^ai---afe^mi 



yCLk |2 



(4) 



where the a indices run over the two components of the 
spinor z. In this representation, the Hamiltonian is a 
simple polynomial, but the coordinates z are constrained 
to be normalized z°(z°)* = 1 and the local phase rotation 
z° — >■ 6*^2;" is redundant. As usual, a convenient choice 
for the spinor z associated to the direction h is 



cos(6l/2) 
sin(6l/2)e*^ 



(5) 



where 6*, (j> are polar coordinates for the direction h. 

Finally, let us define the random ensemble for k- 
P/QSAT in more detail. There are two sources of ran- 
domness in the ensemble: 1. the discrete choice of inter- 
action graph G and 2. the continuous choice of projectors 
n™ associated to each edge. The latter choice is particu- 
larly powerful: generic choices of projectors reduce quan- 
tum satisfiability to a graph, rather than Hamiltonian, 
property [6]. More precisely, for fixed G, the dimension 
Rg — I kei H\ of the satisfying subspace of H is almost 
always minimal with respect to the choice of projectors 
n™. Thus, if G can be frustrated by some choice of pro- 
jectors, it is frustrated for almost all such choices. We 
refer to this as the "geometrization" property. 

Also of particular interest to us is the result that zero 
energy product states for generic QSAT instances exist 
if and only if the factor graph G has 'dimer coverings' 
between its variables and clauses which cover all of the 
clauses [7]. We review this construction in some detail 



below (see Sec. Ill A). Throughout this paper, we use the 
term generic to refer to the continuous choice of projec- 
tors and random to refer to the choice of the graph. 

As for the discrete choice of random interaction graph 
G, we follow previous work and use the Erdos-Renyi en- 
semble with clause density a, in which each of the po- 
tential ( ^ ) edges appears with probability p = aN/ ( j, ) . 
One of the most important geometric features of such in- 
teraction graphs is the existence of a (hyper-)core above 
a critical ahc- The hypercore G" C G is the maximal 
subgraph of G such that every node has degree at least 
two. One can show that all non-core clauses in G — G' 
may be satisfied by local product states, whether or not 
G' is satisfiable, thus all of the nonzero energetics are 
associated with the core. 



A. Low temperature 

Generic PSAT has zero energy satisfying states when 
the interaction graph G possesses 'dimer coverings': 
matchings of spins to clauses such that every clause is 
covered by exactly one dimer and no spin is used more 
than once [7]. Let us review the construction of these 
states. The counting is easiest if we assume that the 
projectors are all of product form 






^t: 



(7) 



where each 2-spinor 0™ is individually normalized. A 
zero energy product state may be associated to each 
dimer covering of the factor graph: if site i is paired with 
clause m by the dimer covering, we use site i to satisfy 
clause m by making it orthogonal to the relevant factor 
(e"'' is the antisymmetric tensor) 



^ab li 



(8) 



This leaves N — M spinors completely unspecified - these 
are zero modes of the system. Meanwhile, any deviation 
of the matched spinor at site i leads to a quadratic en- 
ergy cost for generic choices of the unspecified spinors, so 
that there are M complex harmonic modes, or 2M real 
harmonic modes. Recalling that the ground state energy 
is 0, equipartition produces 



U w TM = TaN. 



(9) 



This counting persists for generic projectors (not of 
product form). To recap [7], let {zi} be a zero energy 
ground state of iJ. Since H{{zi}) = 0, the Zi must satisfy 



III. LOW AND HIGH TEMPERATURE 
EXPANSIONS 

We begin our study of the PSAT model by considering 
the leading order low and high temperature expansions 
for the classical Gibbs ensemble, 



II Dn, 



-PH{{n,}) 



(6) 



These limits provide non-trivial consistency checks for 
the more sophisticated numerical and cavity analyses in 
later sections. As the model is classical, the low temper- 
ature expansion follows from equipartition, even in the 
presence of disorder averaging. However, one must be 
careful to take account of 'zero modes'. A classical zero 
mode corresponds to a degenerate continuous manifold 
of ground states which therefore do not contribute to the 
specific heat. We work this out in Sec. Ill A below. 



The disorder-averaged high temperature expansion is a 
straightforward exercise in combinatorics. We summarize 



the leading order details in Sec. IIIB 



C...a.CV--C. =0 Vm = l...M (10) 

By rotating the local bases, we can assume that the 

ground state Zi = „ J for all sites i, in which case 

these equations simply reduce to ^"..o = 0. Perturbing 
/ 1 



z,; — > I c ) , we find (to leading order) 



j=i 
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SWr, 



(11) 



The dimension of the kernel of this set of M linear equa- 
tions in N unknowns Swi gives the number of (complex) 
zero modes. For factor graphs G with dimer coverings, 
this linear map is generically surjective so the dimen- 
sional of the kernel is A^ — M and we recover the counting 
from the product projector case. 

For graphs that are not product satisfiable, the ground 
state energy is non-zero on the hypercore of the graph G 
and we expect no zero modes associated with the core 
spins. Nonetheless, there may still be zero modes on 
the non-core spins. An equipartition based estimate of 



the low temperature energy density in the non-product 
satisfiable phase looks like: 



[/«Afeo + T(7Vc + (M-M,)) 



(12) 



where eq is the ground state energy per clause (eg > 
0), Nc is the number of spins on the core, Mc is the 
number of clauses on the core. For random interaction 
graphs G, the geometric quantities Nc{a) and Mc{a) are 
known |26| and we have quoted representations for them 
in Appendix [C| 

Of course, this estimate only applies at sufficiently low 
temperatures that all of the non-zero modes indeed look 
harmonic. 



B. High temperature 



Order 


Symbol 


Cumulant 


Value 


Geometry 


1 


^0 


(E")^ 


1x2"^ 


M 


2 


^0 


{(Eor). 


7/9x2-« 


M 


3 


„(3) 
^0 


{(E^r). 


14/15x2"'' 


M 


4 




{(E^'Y). 


5078/4455x2-^2 
-2/243x2^^2 


M 




{(£0)2(i?i)2)^ 


3Mak^ 



TABLE I. The non-vanishing terms in the high temperature 
expansion up to order 4 explicitly for k — 3. The (•)^ indicates 
the symmetrized joint cumulant of the clause energy variables. 
The geometric factor indicates the expected number of such 
terms that contribute to the free energy in the random graph. 
We break out the factors of 2"*" (n being the order) in the 
value as we expect the terms to scale as such. 



The high temperature expansion for the disorder aver- 
aged free energy is straightforward: 
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(13) 



where (•) denotes the infinite temperature average and ~ 
denotes the disorder average. Normalizing the volume of 
the spinor space to 47r, we have 



logZ(O) =7Vlog47r 



(14) 



independent of disorder realization. 

Both the thermal and disorder correlation functions 
which arise in evaluating the high temperature expan- 
sion are with respect to normalized vectors on complex 
spherical spaces. These satisfy a variant of Wick's the- 
orem which we derive in Appendix [B] For example, to 
evaluate the disorder averaged infinite temperature en- 
ergy: 



(^>o = E(^" 



2fe 2 ■" 2 

_ ^^ 

In evaluating the higher order cumulants, two prop- 
erties are especially important. First, as usual in the 



high temperature expansion, disconnected clusters van- 
ish. Thus, the second order cumulant is 



c*^) = Yl (-^™-^"> - (-^'") (^") 



- j^{Hi) - 4 (i73) {H) - 3 (i/2)2 + 12 (i/2) (^Hf _ 6 (ij)'i 



M(4'^+afc2cf^) 



where the constants 
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(15) 



(16) 

(17) 



correspond to the single clause energy fluctuations and 
the overlapping clause energy fluctuations and ak'^ is the 
expected number of clauses intersecting a given clause m 
in one spin. 

Second, the disorder average of any power of a single 
clause energy function produces a number with no ther- 
mal fluctuations. 



{E^iz)y 



2fe(2fc + i)...(2fc-f n-l) 



(18) 



Thus, in any term in the cumulant expansion where all 
copies of a clause lie in the same thermal average, we may 
disorder average first and extract the constant. This im- 
plies that any cumulant in which a particular clause E"^ 
arises only once vanishes, as 'constant' random variables 
always lead to vanishing joint cumulants. For example, 
this tells us immediately that ci ' = 0. 

Using these rules, we find that the only terms that con- 
tribute to the high temperature expansion to 4th order 
are those given in Table |lj We have included the explicit 
evaluation of those terms for k — 3. 

To summarize, the thermodynamic densities to fourth 
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FIG. 2. Finite energy density phase diagram for 3-QSAT ex- 
tracted from various numerical methods. Dark blue is QSAT 
inaccessible region, as estimated by exact diagonalization for 
size A'' = 20. Medium blue region is PSAT inaccessible region, 
as estimated by local quench for size A^ = 100. The light blue 
region indicates where the cavity analysis of PSAT exhibits a 
dynamical instability. More details regarding the numerical 
techniques and averaging are in the main text. 
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IV. FINITE-SIZE NUMERICS 

The general considerations regarding zero energy 
states and zero modes provide a surprisingly detailed 
picture of the satisfiable phase of PSAT, including its 
low energy density of states in both the satisfiable and 
non-satisfiable phases. However, these arguments fail to 
provide any estimate of the ground state energy density 
when it becomes finite, nor do they help to understand 
the landscape of the search problem for an algorithm at- 
tempting to optimize a particular instance. In this sec- 
tion, we turn to several complementary finite-size numer- 
ical studies in order to pin down the PSAT ground state 
energy density and compare it to the QSAT ground state 
energy density. 



A. Greedy local quench 

The most successful numerical line of attack on the 
PSAT ground state energy density turns out also to be 
the simplest: the greedy local quench. The algorithm 
proceeds as follows: Given a PSAT instance H, generate 
a random initial spin configuration 2:^(0) and then itera- 
tively sweep over the sites i in some fixed order. At each 
step, update the local spin to point in the direction which 
optimizes the energy of the neighboring clauses: 



Zi{t + 1) = arg min > 



£;™(z„{z,(i)bea™-.)- (22) 



Here, the integer t keeps track of the number of sweeps in 
the algorithm. The calculation of the local minimizer is 
straightforward. One simply calculates an effective field 
hfn on site i due to each clause m and then orients the 
spinor Zi in the direction of the net field ^^^ h^. As the 
energy function decreases with every step of the quench, 
the algorithm is guaranteed to converge toward a local 
minimum of the energy. We terminate the quench when 
the relative change in the energy SE/E after a full sweep 
across the N sites is less than a fixed threshold (10^^ in 
all data shown). In applying the quench to the random 
ensemble, we have generated several hundred instances of 
size A^ = 100 at varying clause density a, stripped them 
to their hypercores (as all clauses not in the hypercore 
can be satisfied independent of the state of the core) and 
run the quench on each core with ten different initial 
conditions, taken the lowest resulting energy. 

The greedy quench does an excellent job of finding low 
energy product states in the random 3-PSAT ensemble, 
as shown in Figure [2j The quench readily finds states 
with energies equal to zero to machine precision below 
aps — 0.9 ± 0.06 (see width of scatter region in Fig. Is] 
inset). Above this threshold, the found energies begin 
to grow from zero, in agreement with the rigorous ex- 
pectation that the ground state energy is no longer pre- 
cisely zero. We note that the observed error in aps is 
consistent with the fiuctuations in the core density due 
to finite-size N = 100, which we can roughly estimate 

Nc/Mc = 1 ± -^ by noting that N^/N = M^/N = 0.6 



N 



in the thermodynamic limit and assuming that the choice 
of clauses and sites remaining in the core are selected by 
an independent random process from the initial random 
graph with size N. 

The observed energy density for a > aps is both quan- 
titatively very small - even at a = 8, the found energies 
of 0.05 per clause are significantly below the naive high 
density estimate of 1/8 — 0.125 per clause - and the 
critical turn on for a near aps appears very soft (see in- 
set, Fig. [2|. The simplest rigorous upper bound for the 
large a energy density is to satisfy the first apgN of the 
clauses and then upper bound the energy of the remain- 
ing clauses by 1/8. This leads to E/M < (1 - aps/a)/8. 

A better estimate at large clause density can be made 
by noting that every spin i is subject to a mean exchange 
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FIG. 3. Product state ground state energy found by greedy 
local quench near critical region a^s ~ 0.92 at finite size 
A'^ = 100. Error bars indicate instance-to-instance fluctu- 



ations over Nsan 



3 — 400 instances per density alpha. 



(inset) Semilogarithmic scatter plot of same data. 



FIG. 4. Typical energy traces as a function of time for 10 
independent sweeps in a typical instance at each of a = 
0.88, 0.90, 1.0 (iV = 100). The curve in the left panel continue 
to fall exponentially (linearly on the log-scale) to machine 
precision (10~^^). 



field due to the di clauses attached to it, where the degree 
di is a Poisson distributed random variable with mean 
ka. With respect to an infinite temperature spin state, 
this exchange field is a sum of di independent random 
variables so that it has 0{\/d'i) ^ 0{\/ka) fiuctuations. 
Thus, by aligning each local spin with its local exchange 
field, we expect to gain of order 0{-\/ka) energy per spin 
relative to the infinite temperature state. This produces 
a high density estimate of E/M = 1/2'= - 0{l/y/a). We 
believe this essential scaling at large a is robust although 
the coefficient hiding in the 0-notation is not so easy to 
produce - the variance of the exchange field in an in- 
finite temperature state is straightforward to calculate, 
but the process of tipping spins toward their local ex- 
change field modifies the field on neighbors in a nontrivial 
manner. In any event, we note that a fit (not shown) to 
this quadratic form at high clause density a > 2 produces 
good agreement with an extracted coefficient of roughly 
0.18-0.21/VS (at A: = 3). 

The dynamics of the local quench show mild evidence 
for the presence of metastable minima in the classical 
PSAT energy landscape above the transition region, as 
shown in the typical time traces shown in Fig. l4J For 
densities below aps = 0.91, the typical traces decay con- 
tinuously to very low values while at higher density, they 
exhibit a larger variance of found energies for different 
initial conditions and step-like metastable behavior in the 
individual trajectories. 



B. Simulated annealing 

As a second approach to finding low energy product 
states, we implemented a standard simulated annealing 
protocol to cool the PSAT Hamiltonian to low temper- 
ature on randomly generated instances. At each step of 



this procedure, a randomly sampled spin is perturbed by 
a small random step and the move is accepted following 
a Metropolis rule at temperature T . The temperature 
decays exponentially with a decay constant of 1.003 per 
10'* Monte Carlo sweeps, beginning with temperature 1 
and finishing at T = 10~*. 

The resulting disorder averaged energies A^ := 10 — 80 
qubits show extremely limited finite size scaling (not 
shown) so that we have included only the largest size 
data {N = 80) in Fig. [2] We note that simulated an- 
nealing fails to find energies as low as the local quench 
protocol but agrees well with the low temperature ex- 
trapolation of the cavity results (see Sec. W\. While it 
is possible that the annealing results would improve by 
further optimizing the cooling schedule, we suspect that 
the timescale for this cooling may become very long in 
the UNSAT regime. This is indicated by the transition in 
the behavior of the step size dependence of the simulated 
annealing results as the clause density is tuned from the 
SAT to UNSAT phase, see Fig. [5] 



C. Exact diagonalization 

Finally, we include exact diagonalization results for 
systems of up to A^ = 20 qubits in Fig. [2J These pro- 
vide a numerical baseline for the QSAT ground state en- 
ergy density in the UNSAT regime (striped symbols) and 
demonstrate its separation from the best found product 
states by local quenching. We note that the finite size 
scaling of these results looks quite converged for /c = 3 
(not shown), as seen in previous studies as well [7], but 
that for larger k there are still strong finite size effects 
for the achievable system sizes, as the graphs are quite 
unlike the thermodynamic random graph at these sizes. 
This is indeed the central reason that we limit the nu- 
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FIG. 5. Step size dependence of the best found PS AT energy 
in the simulated annealing algorithm. The step size controls 
the maximum change of a single component of the wave func- 
tion per update trial. In the PRODSAT phase reducing the 
step size leads to a monotonically decreasing energy. In the 
UNPRODSAT phase the behavior is irregular. 



merical investigations of this paper to the k — 3 case. 



V. CAVITY APPROACH 

In this section, we look at the PSAT model using so- 
called cavity techniques [11]. Our motivations for do- 
ing so are twofold. Firstly, while the previous finite size 
numerical work obtained low ground state energy esti- 
mates, these were limited to relatively small system sizes 
of at most N = 100 variables. With the cavity approach 
we can explore the thermodynamic limit where iV — > cx) 
directly. Secondly, we also observed regions of the pa- 
rameter space where the finite size algorithms began to 
slow down. In the local field quench, this slow-down coin- 
cided with what appeared to be metastable states, while 
in the simulated annealing results, the algorithm reached 
a barrier in finding lower energy ground states (see Fig- 
ures [5] and [4]). However, these observations in a finite size 
system do not indicate whether these metastable states 
refiect true transitions in the thermodynamic structure 
of the configuration space. 

The cavity approach allows us to detect the appear- 



ance of such states as A^ — >■ oo. This is indicated by 
the appearance of a dynamical instability where the sim- 
plest (replica-symmetric) cavity ansatz breaks down [T^] . 
From a algorithmic perspective, it is unclear whether 
such a boundary is truly relevant, as indicated by the 
robust results of the local field quench. However, in the 
spin glass lore, it is believed that such barriers ultimately 
limit local search procedures for discrete classical systems 
and we believe the discovery of similar instabilities in the 
continuous degree of freedom PSAT system represents an 
interesting avenue for future detailed exploration. 

In the next section, wc introduce the cavity approach 
appropriate to the PSAT problem and provide the ther- 
modynamic quantities that can be obtained from it. Af- 
ter this, we discuss how the replica symmetry broken 
phase can be detected in the cavity approach. 



A. Cavity method 

The cavity approach to the classical PSAT model con- 
sists of studying the cavity distributions Pi^mini). This 
is the marginal classical probability distributions for the 
unit vector fii at site i, in the absence of the clause m. 
The so-called belief propagation equations follow from 
assuming that the cavity distributions for the neighbors 
of the clause m are independent when m is removed. Un- 
der this assumption, the cavity distributions as shown in 
Fig [6] obey the equations 



j^drn—i 



Pi^m{ni) = ^ Yl Qn^^i'^i) 



(23) 

(24) 



n^di — m 



Here, Q{n) is an intermediate cavity function that is use- 
ful primarily as a bookkeeping device. For a /c-local in- 
teraction graph with M clauses, this provides 2AIk func- 
tional equations for 2Mk unknown cavity distributions 
P{n) and Q{h). 

In the thermodynamic limit, we expect that the cavity 
distributions themselves are i.i.d random variables with 
a self-averaging distribution P [P(-)] and P [Q(-)]- Under 
this assumption, we obtain self-consistent equations for 
the functional distributions P [P(-)] and P [Q{-)]: 



'[F(-)]=Ed,Q„<5 



Pi-)-^f[Qn{-) 






(25) 
(26) 



These are the so-called replica-symmetric cavity equa- 
tions for the PSAT model. Here, E stands for an expec- 
tation value with respect to the subscripted quantities, 
such as d, the onward degree of a spin, and E™, the choice 
of projector on the clause. The normalization factors Z 
implicitly depend on these random variables. 




FIG. 6. The definition of the cavity distributions on a locally 
tree-like interaction graph. Squares indicate clauses, labeled 
by letters, and spheres the spinor variables, labeled by inte- 
gers. 



In general, the cavity equations (25) can not be solved 



exactly for the functional distributions P [P(n)], P [Q(7i)]. 
While the relevant cavity equations are often analytically 
intractable even in models with discrete local degrees of 
freedom, such as in classical fc-SAT, in PSAT, the cavity 
messages P{'h) and Q{h) are general normalized proba- 



bility distributions on the sphere, and the cavity distri- 
butions P [•] are distributions on that functional space. In 
the presence of sufficient symmetry, as, for example in ro- 
tationally invariant models [TT, one can make progress 
in the symmetric phase analytically, but PSAT has no 
such rotational symmetry and we must resort to numer- 
ical methods throughout the phase diagram. 

Thus, we follow established practice[TTJ [?7H^ in 
studying cavity equations by representing the distribu- 
tions, P [P{-)] and P [Q(-)]; by a large population of sam- 
pled messages P(-) and Q{-)- These probability distribu- 
tions on the sphere in turn are represented by discretized 
probability vectors Pi, with a finite number i — 1 ■ ■ ■ iVdisc 
of patches roughly equally sampled over the surface of the 
sphere. By varying the discretization, we can check for 
convergence of thermodynamic observables to the contin- 
uous degree of freedom iVdisc ^^ oo limit. 



To solve the cavity equations (251, we employ popula- 
tion dynamics. We initialize the population of messages 
P(-) and Q{-) randomly and then proceed iteratively to 
use the belief propagation equations ( p3| ) to generate the 
new messages from the existing population. As the pop- 
ulation size Npop and total number of time steps T ap- 
proach infinity, we expect the resulting population to 
converge to a fixed point of the cavity equations jllj . 
Of course, in practice, we take a fixed population size of 
Npop — 10^ and iterate until we observe statistical con- 
vergence of various moments. 

Finally, with the converged messages in hand, one can 
estimate the various disorder averaged thermodynamic 
quantities such as the internal energy U and the Bethe 
free energy, F. 



For example, the internal energy per clause follows directly from the expected value of the energy of a randomly 
added clause in the presence of k sampled cavity messages. 



U/M ^Esa^P 



I [ Dn,Dn2DhsE-{n,,n2,n3)e-^'''^'''^^-'^''>P,^a{ni)P2^a(n2)P3^a{n3) 



(27) 



The expression for the Bethe free energy per variable / — F/N is somewhat more complicated [H] , as it follows from 
balancing the change in free energy of various graph surgery operations. 



/(a,^) = fv{a,/3) + afa{a,P) ~akfe{a,l3) 



where, 



(28) 



-/3A=E,,Qlog 

-/3/a=E£a^plOg 

-/3U=Ep^Qlog 



DfiQi^i{h) . . .Qe^i{h) 

DniDn2Dn3e-^''''(-'''^^-'''-''>Pi^a{ni)P2^a{n2)P3^a{n3) 
DhP{n)Q{h) 



(29) 



These represent the change in the free energy from the 
addition of a single clause which has contributions com- 



r 



ing from the variables added (/t,), the clause itself (/a) 
and each edge connecting variable to edge (/g). Of 
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FIG. 7. Cavity estimate of the internal energy of PSAT as 
function of temperature with Ndiac = 24 states per Bloch 
sphere. Error bars reflect samphng error in calculating E 
from the converged population. Solid lines indicate the high 
temperature expansion of the energy to 4th order. Solid (hol- 
low) markers indicate data in the replica symmetric (broken) 
phase, (inset) Zoomed in at low temperature of same data, in- 
cluding results for both Ndiac = 24 states (■) and Ndiac = 12 
states (•). Solid (dashed) lines indicate T-linear fits to the 
Ndiac = 24 {Ndrsc = 12) data. 



course, with the internal energy IJ and free energy F 
in hand, it is straightforward to extract the entropy S 
from the thermodynamic relation F — U — TS. 



B. Cavity thermodynamic results 

We plot the energy density in Fig. [7] and free energy 
density in Fig. [8] at various temperatures, clause densi- 
ties and discretizations. The simplest results are at high 
temperature where the discretization error becomes small 
(not shown) and we expect detailed agreement with the 
high temperature expansion. Indeed, the expansion to 
4th order, given by equations dloj 
well for /3 < 



(21 ), works extremely 
as can be seen from the (no fit) 



5-10 

agreement between the data and the solid lines in both 
the energy and free energy plots. We view this agreement 
as confirmation of both the high temperature expansion 
and the validity of the numerical cavity analysis. 

At low temperature, we can fit the cavity results to the 
equipartition analysis of Sec. |III A| In particular, we ex- 
pect the energy to exhibit T-linear behavior with a zero- 
temperature intercept giving the ground state energy of 
the model and a slope corresponding to the specific heat. 
We have converged cavity populations down to temper- 
atures of T = 1/500 as shown in the inset to Fig. [t] 
The linear extrapolations provide an estimated ground 
state energy of roughly 10""* for a = 0.25, 0.75, 1.00 while 
for a = 1.25,2,3,4, the ground state energy is of order 
10~^. These intercepts can be seen in Figure [2] We note 
that the statistical error bars in these estimates are quite 



FIG. 8. The free energy, F/N as a function of a, T for qubits 
that have been discretized into 24 states. Point labeling is the 
same as in Figure [7] and likewise, hollow points indicate data 
collected ignoring the dynamical replica instability, (inset) 
As a function of /3, we show the fit to the high temperature 
(small /3) approximation for the free energy (cf. equations 
(19|- pll). 



large relative to the data found in the local quench nu- 
merics and, although they in principle exhibit no finite 
size effects, they do contain discretization error, which 
we expect to systematically raise the estimated energies 
and suppress the low temperature specific heat, as can 
be seen in the inset to Fig. [7] Also, as we will discuss 
in more detail below, the low temperature-high density 
cavity calculation exhibits a dynamical instability below 
which we should hesitate to believe the quantitative re- 
sults. Nonetheless, the ground state energies predicted 
by the cavity calculation at Ndisc = 24 are consistent 
with those found by simulated annealing and, at low den- 
sity, with the zero expected for a < aps — 0.91. 



C. Dynamic instability 

In the context of discrete classical systems, the cav- 
ity equations ( 25 1 exhibit a well studied collection of 



instabilities that are believed to indicate the prolifer- 
ation of solutions to the underlying belief propagation 
(BP) equations (23) on large instances [TTJI3D]. In turn, 
each of the many, extensively distinct solutions to the 
BP equations has been argued to correspond to a dis- 
tinct thermodynamic pure state in the Gibbs ensemble. 
In the replica treatment of spin glass models, these in- 
stabilities are associated with replica symmetry break- 
ing, and this language has been largely imported into 
the cavity formalism. From an algorithmic point of 
view, replica symmetry breaking indicates the appear- 
ance of metastable states that are 'thermodynamically 
separated'. Such metastable configurations are believed 
to trap local search and annealing algorithms, frustrat- 
ing the approach to lower temperatures - a behavior that 
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we observed in our finite-size numerics in the high clause 
density regime. 

Fortunately, while the full treatment of replica symme- 
try breaking requires the solution of a yet more compli- 
cated hierarchy of equations than (251, it is possible to 



detect the presence of replica symmetry breaking within 
the replica symmetric treatment [311 132) . To perform 
such a study, one initializes the population dynamics 
algorithm with two populations that differ by a small 
random displacement, 6V, and then evolves each popu- 
lation with an identical sequence of random moves. In 
the replica symmetric phase, these two initial conditions 
result in identical distributions for Pi^m and Qm^i bX 
late times. In the presence of replica symmetry breaking, 
the two distributions diverge. The Lyapunov exponent 
associated with this process (inset. Fig. |9]) provides a ro- 
bust indication of the underlying instability: where the 
exponent goes from negative to positive, we identify the 
dynamic instability phase boundary in the (a,r) plane 
(Fig. |9| . As can be seen best in Fig. [2] the high energy, 
low clause density regime of PSAT exhibits full replica 
symmetry while the dynamical instability sets in along 
a boundary which encloses the low energy density, high 
clause density regime, frustrating the approach to the 
ground state in the non-PRODSAT phase. 

As our numerical approach relies on discretization of 
the continuous spins, we have attempted to check the 
convergence of the instability boundary to a finite limit 
as Ndisc becomes large. Indeed, for any discretization 
we find that there is an ac{T) where there is a dy- 
namical glass transition. As the discretization is in- 
creased, the instability line moves to higher clause density 
but, for the three discretizations probed, we believe this 
finite-discretization shift converges to a finite tempera- 
ture boundary for a > 0.9 (see Fig. [9|. This is especially 
striking for, as T — > 0, we find that this transition ap- 
parently coincides with the zero temperature PRODSAT 
transition. [14j We note that this is different from the sit- 
uation in XORSAT where the dynamical glass transition 
at zero temperature coincides with the appearance of the 
hypercore at Uhc ^ 0.81 rather than the SAT-UNSAT 
transition at a — 0.91 1^. 



VI. CONCLUDING REMARKS 

As noted at the outset this paper is fundamentally 
concerned with the approximation problem for randomly 
chosen QSAT instances. Specifically, by combining sev- 
eral methods — cavity analysis, simulated annealing and 
a greedy algorithm — we give an estimate of the limits to 
which product states can be used to approximate solu- 
tions to 3-QSAT. Among these methods, the cavity anal- 
ysis shows a dynamical glass transition before the limits 
of product state approximation are reached but this does 
not appear to have a marked influence on the perfor- 
mance of our best algorithm, the greedy search - at least 
for sizes up to TV = 100. As such it is not clear whether 
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FIG. 9. (inset) The sum-of-squares difference between two 
populations as a function of time in the population dynamics 
for Ndiac = 36 states and /3 = 100. From bottom to top in the 
inset, we show data for a — 1.26, 1.28 and 1.29. We observed 
a transition in a between when the populations converge or 
diverge over time. This boundary is plotted in the main fig- 
ure, (main) The location of the dynamical instability phase 
boundary extracted from the Lyapunov exponents of the pop- 
ulation dynamics. As the discretization of the Bloch sphere 
increases, the dynamical instability appears to converge to a 
finite temperature phase boundary. As T — >■ 0, the instabil- 
ity line appears to terminate in the PRODSAT transition at 
a ~ 0.91. 



the search for product states exhibits an easy-hard tran- 
sition before the limits of the product state region are 
reached. Along the way, we have provided various addi- 
tional statistical mechanical results on the classical PSAT 
Hamiltonian which may be of of interest to aficionados 
of disordered systems. 

A natural next step is to study classical wavefunctions 
that build in some entanglement. It may be possible to 
analyze bounded-depth local quantum circuits matched 
to the random graph structure, similar to those studied 
in [20]. It is unclear how the non-commuting nature of 
generic QSAT Hamiltonians affects the quality of this 
approach and whether one needs to move to more general 
entanglement structures. We expect such considerations 
and the set of conjectures embodied in our opening figure 
(Fig. 1) to provide stimulus for future work. 
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Appendix A: Exact singletons 

1. Singleton clause 

The partition function for a single k ~ \ clause at- 
tached to a spin is given by: 



Z ^ / d{cos9)d(j)e 



'/3|0a2°|' 



(Al) 



Choosing the coordinate system to align with the state 
16), this reduces to 



Z 



d(cos0)#e-^-^ W2) 
27r /"d(cos0)e-^(i-^°^^)/2 



27re 



^/3/2 



./9«/2 



13/2 
47r(l - e-^)/l3 



(A2) 
(A3) 

(A4) 
(A5) 



The free energy is 

F = -- 

I 

The entropy is 
OF 



log 4tt + log 






S 



dT ^ d(3 
= log An + log h 1 — 



3/3-1 



The energy is 



^^^ + ^/^4-e^- 



(A7) 
(A8) 

(A9) 



These expressions recover the correct limits (high tem- 
perature): 
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C/^-(l 
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13'^ l + /3/2 + ^73!' 2 



S^\ogATT + Oi/3^) 
and low temperature: 



/3/6 + 0(/32) 

(AlO) 
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(A12) 

5'-log47r//3 + l + C'(e-'^) (A13) 

To check the high temperature expansion, we have 

'l-u\ I 

,2 ' 



{H), = {\<i,aZr)^ = l^- 
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2-^2 



(A14) 
(A15) 
(A16) 



Here, we note that the infinite temperature measure on 
u is simply uniform on [—1, 1]. 



2. Single product clause 

For a product projector, the integrals for the partition 
function may be simplified by rotating the local bases to 
match those of the projector. Thus, 

Z = f d{cos ei)d(j)i ■ ■ ■ d{cos 9k)d(j)k g-'^l'^-^i I'-l-^^^^l" 

(A17) 

== (27r)'= f dui--- dukc'^'-^-'-^ (A18) 

For fc = 2, this can be evaluated to the useless form, 

Z=^[7 + r[0,/3]+log/3] (A19) 



Appendix B: Uniform measure on 



(A6) The measure on the projectors is the uniform Haar 



measure on the space ' 



We can represent this mea- 



sure on an n-component complex vector 
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With the partition function in hand we can evaluate 
the disorder averaged correlators. The two-point corre- 
lator is 
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The higher order correlators satisfy a sHghtly modified 
Wick's theorem: 



'^.n - M^L±M^ (B2) 



K4>2 • • • 0fc0fe+10/c+2 ■ ■ ■4>2k 



n{n + 1) 

Err c 

pairi7igs 

n{n -\-l) ■ ■ • (n -\- k — 1) 
(B3) 



These averages are equaUy useful for evaluating high 
temperature correlators of the two-component spinors Za- 



where the average degree A* corresponds to the largest 
solution of the equation: 



e-^' - 1 



ka 



= 0. 



(C2) 



The total number of nodes in the hypercore is 



Appendix C: Hypercore data 

We summarize a few results derived in |26| regarding 
the hypercore of random factor graphs G. The hypercore 
is the maximal subgraph of G on which all nodes have 
degree at least 2. The quoted results follow from ana- 
lyzing the leaf removal algorithm applied to the random 
graph. In the notation of 26 , the clause density a is 7 
and the number of spins per clause k is p. 

The degree distribution on the hypercore is Poissonian 
for degrees d > 2: 



for d = 0, 1 



for d>2 



(CI) 



Nc{a) = N^P^{d) = N \l - {I + X')e 



,-A* 



(C3) 



d>2 



and the total number of edges (clauses) 



M,(a)=iV-(l~e-^*). 



(C4) 



For large a, the hypercore takes over most of the graph 
up to exponentially small corrections: A* ~ ka and N^. — 
N{1 -{1 + ka)e-''°'), M^ = Na(l - e"*^"). 
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